

%%
fs     = Head.SampRate;
dur    = 30;

F      = freq_GetFreqBins(fs*dur,fs);
F      = F(1:16:end);
fidx   = util_GetIndex(F,10):util_GetIndex(F,16);
freq   = F(fidx);

chan   = Head.DataChan;
n_chan = length(chan);
n_time = Head.Stage.Duration * fs;
n_epch = length(Head.Stage.Time);

stime  = Head.StartDate + Head.StartTime;
sidx   = util_GetOffIndex(stime,fs,Head.Stage.Time(1));
eidx   = sidx+(n_time * n_epch)-1;

n2     = util_GetIndex(Head.Stage.Label, {'N2'});
epch   = ~isnan(Head.Stage.Series(:,n2));

% wavelet
deley   = fs * 0.4;
spindle = [];
c       = util_GetIndex(Head.ChanLabel, {'C3'});
e       = 1;
thr     = mean(abs(Data(c,sidx:eidx)),2) * 4.5;


%%
fs = Head.SampRate;
Wn = [10 16];
win = 'hann';
mul = 5;
c = util_GetIndex(Head.Stage.Label, {'N2'});

fData = FirFilter(Data(c,:),fs,Wn,win,mul);
C = abs(hilbert(fData));
Z = zscore(C);
peaks = getpeaks(Z,1);
tmp = find(Z(peaks) > 2.328); % one-tail P = 0.01
candi = peaks(tmp);



for e = 1 : length(epch)
    if epch(e)
        etime = Head.Stage.Time(e);
        sidx = util_GetOffIndex(stime,fs,etime);
        eidx = sidx+6000-1;
        [~,idx] = find(bitand(candi >= sidx, candi <= eidx));
        spindle = [spindle candi(idx)];
    end
end




%%
s = 1;
spdl = abs(lgt_cwt(Data(c,spindle(s)-100:spindle(s)+100)',fs,freq,9,'gab','same'));
figure;imagesc([],freq,spdl);

t = util_GetOffTime(Head.StartDate+Head.StartTime,Head.SampRate,spindle);
out = datestr(t,'HH:MM:SS.FFF');
xlswrite('spindle.xlsx',out);












































%%

for i = sidx : n_time : eidx
    
    if epch(e)
        seg        = Data(c,i:i+n_time-1);
        E          = mean(abs(lgt_cwt(seg',fs,freq,9,'gab','same')),1);
        [~,candidate] = find(E > thr);
        
        for d = candidate
            if all(E(d:d+deley) > thr)
                spindle = [spindle i+candi-1];
            else
            end
        end
        
        
        for j = deley/2+1 : n_time-(deley/2)
            
        end
        
        
        
%         
%         peaks      = getpeaks(E, 1);
%         candidate  = peaks(E(peaks) > thr);
%         
%         for candi = candidate
%             lidx = cand-deley;
%             ridx = candi-deley;
%             
%             if lidx < 1
%                 lidx = 1;
%             end
%             if ridx > n_time
%                 ridx = n_time;
%             end
%             
%             while ridx-lidx >= deley
%                 if all(E(lidx:ridx) > thr)
%                     spindle = [spindle i+candi-1];
%                 else
%                     
%                 end
%             end
%             
%             if lidx > 0 && ridx <= n_time
%                 if all(E(lidx:ridx) > thr)
%                     
%                 end
%             end
%         end
           
%         [~,idx] = max(E);
%         lidx = idx-deley;
%         ridx = idx+deley;
%         if (lidx > 0) && (ridx <= n_time)
%             if all(E(lidx:ridx) > thr)
%                 spindle = [spindle i+idx-1];
%             end
%         end
         
    end
    e = e + 1;
end

%%
s = 2;
spdl = abs(lgt_cwt(Data(c,spindle(s)-100:spindle(s)+100)',fs,freq,9,'gab','same'));
figure;imagesc([],freq,spdl);

t = util_GetOffTime(Head.StartDate+Head.StartTime,Head.SampRate,spindle(s));
datestr(t,'HH:MM:SS.FFF')


% Head.Event.Label = {'Spindle'};
% Head.Event.Color = {'R'};
% Head.Event.Time = util_GetOffTime(Head.StartDate+Head.StartTime,Head.SampRate,spindle);
% Head.Event.Series = ones(1,length(spindle));
% Head.Event.Latency = [];